Marxan in R
ESM 270 Tutorial
Cameryn Brock
6/27/2020
Introduction
Designing reserve networks requires considering both ecological and socioeconomic factors. Marxan is a commonly used spatial planning tool that quantifies these factors to generate efficient reserve networks. Here, we will use species and parcel data from a 2014-15 Bren Group Project in the Morro Bay Watershed to examine the optimate regional reserve network.
Priotitizr
Prioritizr is an R package that, like Marxan, can guide systematic reserve design by solving conservation prioritization problems. In this tutorial, we will be using the marxan_problem() function which can read Marxan input data. More information on the prioritizr package can be found at https://mran.microsoft.com/snapshot/2018-03-04/web/packages/prioritizr/vignettes/prioritizr_basics.html.
A really useful tutorial in using prioritizr, specifically in using summed solutions, can be found here: https://cran.r-project.org/web/packages/prioritizr/vignettes/tasmania.html. This resource helped frame this tutorial.
Gurobi
For the solve() function, you must have a solver installed. For the the process with summing and multiple runs process, this must be gurobi. This requires following the instructions below.
- Download gurobi (https://www.gurobi.com/downloads/gurobi-optimizer-eula/)
- Request an academic license (https://www.gurobi.com/downloads/end-user-license-agreement-academic/)
- Verify license in terminal by copying and pasting information in your acount > your licenses
- Must install r package with
install.packages('/Library/gurobi902/mac64/R/gurobi_9.0-2_R_3.6.1.tgz', repos=NULL)(if on a mac and downloaded the most recent version of gurobi, don’t download newest version of R! (4.0.0; at least as of 5/2020 this was not functional)). - must install slam package with `install.packages(“slam”)
- add gurobi and slam libraries
Data
We will be using a case study (and the data they prepared) from a past group project on the Morro Bay National Estuary Program.
Species data
- MorroBay_species_polyons: shapefiles of species locations
- MorroBay_species_pts: point data of species locations
- MorroBay_habitats: shapefile of location of key habitats
- morro-bay-spec.csv: information on species/habitat ID (the conservation feature), target value for that species (these values are 30% of the # of parcels containing that species or 10% for habitats), and the species name. The ‘spf’ column is the ‘species penalty factor’ and is a multiplier that determines the penalty if the conservation objective is not met.
- spec-name-status.csv: information on species status.
Parcel data
- morro-bay-pu.csv: includes information on the planning unit ID, the cost of that planning unit, and each parcel’s status (value 0 = available for prioritization, value 3 = locking in, 3 = locked out)
- morro-bay-puvspr.csv: includes information on the amount (presence/absence) of each species in each planning unit (pu)
Getting started
Attach packages
library(here)
library(janitor)
library(prioritizr) # marxan
library(sf) # spatial features
library(slam) # to use gurobi
library(gurobi) # solver
library(ggmap) # basemaps
library(patchwork)
library(tidyverse)If you do not have any of these packages installed, run the following in the console: install.packages("package name")
Note: this will not work for installing gurobi. Instructions are in the introduction section.
Read in data
All .csv files in this tutorial are unaltered from the .xlsx files provided in the R drive. They have only been renamed and saved as csvs. The “data” folder is identical to the MorroBay_data folder provided in the R drive, it has only been renamed.
Species and pu information
spec <- read_csv(here("data", "morro-bay-spec.csv")) %>%
head(140) %>% # mine reads in extra blank rows - skip if yours does not
rename("amount" = "target") # target is called "prop" (relative) or "amount" (absolute)
pu <- read_csv(here("data", "morro-bay-pu.csv")) %>%
select(1:3) # mine reads in an extra blank column - skip if yours does not
puvsp <- read_csv(here("data", "morro-bay-puvspr.csv")) %>%
select(1:3) %>%
head(11849)
status <- read_csv(here("data", "spec-name-status.csv")) %>%
select(1:3) %>%
head(140)Polygons
parcels <- read_sf(dsn = here("data"), layer = "MorroBay_parcels") %>%
clean_names()
ggplot(data = parcels) +
geom_sf()
Analyze all species
Run Marxan
The function marxan_problem() in the prioritizr package gives us the “canned” approach that works for our purposes. If more customizations are desired, feel free to explore the problem() function instead.
Arguments for marxan_problem (for more information, see ?marxan_problem):
- x: planning units to use in reserve design and their cost (
x = pu) - spec: conservation features. Must contain columns for id, name, and prop/amount (
spec = spec) - puvspr: the amount of each feature in each planning unit. Must contain columns for pu, species, and amount. (
puvspr = puvsp) - bound: boundary data (we will not use this;
bound = NULL)
Also note that in solver, we will use 100 runs (number_solutions = 100)
marxan_1 <- marxan_problem(x = pu,
spec = spec,
puvspr = puvsp,
bound = NULL,
blm = 0)
marxan_1_problem <- marxan_1 %>%
add_gurobi_solver(gap = 0.15) %>%
add_pool_portfolio(method = 2, number_solutions = 100)
# gap = optimality gap. 0.15 gives us results closest to that of using marxan software. Learn more under ?add_gurobi_solver
# method 2 finds a specified number of solutions that are nearest to optimatelity. Learn more under ?add_pool_portfolio
marxan_1_soln <- solve(marxan_1_problem)## Gurobi Optimizer version 9.0.2 build v9.0.2rc0 (mac64)
## Optimize a model with 140 rows, 670 columns and 11849 nonzeros
## Model fingerprint: 0xd2c5bccb
## Variable types: 0 continuous, 670 integer (670 binary)
## Coefficient statistics:
## Matrix range [1e+00, 1e+00]
## Objective range [2e+00, 5e+06]
## Bounds range [1e+00, 1e+00]
## RHS range [3e-01, 2e+02]
## Found heuristic solution: objective 4.238988e+07
## Presolve removed 116 rows and 205 columns
## Presolve time: 0.01s
## Presolved: 24 rows, 465 columns, 3066 nonzeros
## Variable types: 0 continuous, 465 integer (465 binary)
## Presolve removed 1 rows and 0 columns
## Presolved: 23 rows, 465 columns, 3058 nonzeros
##
##
## Root relaxation: objective 2.840941e+07, 14 iterations, 0.00 seconds
##
## Nodes | Current Node | Objective Bounds | Work
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## * 0 0 0 2.840941e+07 2.8409e+07 0.00% - 0s
##
## Optimal solution found at node 0 - now completing solution pool...
##
## Nodes | Current Node | Pool Obj. Bounds | Work
## | | Worst |
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## 0 0 - 0 - 2.8409e+07 - - 0s
## 0 0 - 0 - 2.8409e+07 - - 0s
## 0 2 - 0 - 2.8409e+07 - - 0s
##
## Explored 1958 nodes (383 simplex iterations) in 0.12 seconds
## Thread count was 1 (of 4 available processors)
##
## Solution count 100: 2.84094e+07 2.84097e+07 2.842e+07 ... 3.34128e+07
##
## Optimal solution found (tolerance 1.50e-01)
## Best objective 2.840941400000e+07, best bound 2.840941400000e+07, gap 0.0000%Sum solutions
marxan_1_ssoln <- marxan_1_soln %>%
mutate(sum = rowSums(.[6:105])) %>%
select(id, cost, status, locked_in, locked_out, sum)
hist(marxan_1_ssoln$sum,
main = "Selection Frequencies",
xlab = "Number of runs unit was selected")
Create map
Join polygons with output
parcels_marxan_1 <- inner_join(parcels, marxan_1_ssoln, by = "id")
ggplot(data = parcels_marxan_1) +
geom_sf(aes(fill = sum),
color = "white",
size = 0.05) +
scale_fill_gradient(low = "slategray2",
high = "navy") +
labs(fill = "Summed \nSolution") +
theme_minimal()
Add basemap
This step is optional, but will make your map look better and is good to learn for further spatial analyses in R. Unfortunately, the options for basemaps in R are limited. Here, I have used the ggmaps package because it is the most accessible, however, it is still not very accessible compared to other R packages.
Using ggmap requires getting an API key, which is a bit of a hassle, but worth it. Instructions are here: https://cran.r-project.org/web/packages/ggmap/readme/README.html
For a ggmap cheatsheet, including the different basemaps in ggmaps, see: https://www.nceas.ucsb.edu/sites/default/files/2020-04/ggmapCheatsheet.pdf
morrobay <- get_map(location = c(lon = -120.7665, lat = 35.335),
zoom = 12,
maptype = "terrain-background", # - background means no references, omit if want references
source = "google")
all_spec <-
ggmap(morrobay) +
geom_sf(data = parcels_marxan_1,
aes(fill = sum),
color = "white",
size = 0.1,
alpha = 0.85,
inherit.aes = FALSE) +
coord_sf(crs = st_crs(4326)) +
scale_fill_gradient(low = "slategray2",
high = "navy") +
labs(title = "All Species",
fill = "Summed \nSolution",
x = NULL,
y = NULL) +
theme_minimal()
all_spec
Analyze endangered species
Create dataframes
Here we create spec and puvsp data frames similar to those we had for all species, but containing only endangered or threatened species.
end_status <- status %>%
mutate("endangered" = case_when(
str_detect(status, pattern = "endangered") == TRUE ~ "yes",
str_detect(status, pattern = "threatened") == TRUE ~ "yes",
T ~ "no")) %>%
filter(endangered == "yes")
end_spec <- merge(end_status, spec, by = "id") %>%
select(id, amount, spf, name.x) %>%
rename("name" = "name.x")
puvsp_id <- puvsp %>%
rename("id" = "species")
end_puvsp <- merge(end_status, puvsp_id, by = "id") %>%
rename("species" = "id")Run Marxan
Note the following parameters (same as for all species):
- No boundary length (
blm = 0) - Number of runs = 100 (
number_solutions = 100)
marxan_end <- marxan_problem(x = pu,
spec = end_spec,
puvspr = end_puvsp,
bound = NULL,
blm = 0)
marxan_end_problem <- marxan_end %>%
add_gurobi_solver(gap = 0.15) %>%
add_pool_portfolio(method = 2, number_solutions = 100)
# see gap meaning under ?add_gurobi_solver: 0.15 gives us results closest to that of using marxan software
# see method meaning under ?add_pool_portfolio
marxan_end_soln <- solve(marxan_end_problem)## Gurobi Optimizer version 9.0.2 build v9.0.2rc0 (mac64)
## Optimize a model with 18 rows, 670 columns and 970 nonzeros
## Model fingerprint: 0x6e88f508
## Variable types: 0 continuous, 670 integer (670 binary)
## Coefficient statistics:
## Matrix range [1e+00, 1e+00]
## Objective range [2e+00, 5e+06]
## Bounds range [1e+00, 1e+00]
## RHS range [9e-01, 1e+02]
## Found heuristic solution: objective 4.174920e+07
## Presolve removed 14 rows and 203 columns
## Presolve time: 0.00s
## Presolved: 4 rows, 467 columns, 495 nonzeros
## Variable types: 0 continuous, 467 integer (467 binary)
## Found heuristic solution: objective 3.534249e+07
## Presolved: 4 rows, 467 columns, 495 nonzeros
##
##
## Root relaxation: objective 2.677119e+07, 1 iterations, 0.00 seconds
##
## Nodes | Current Node | Objective Bounds | Work
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## * 0 0 0 2.677119e+07 2.6771e+07 0.00% - 0s
##
## Optimal solution found at node 0 - now completing solution pool...
##
## Nodes | Current Node | Pool Obj. Bounds | Work
## | | Worst |
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## 0 0 - 0 - 2.6771e+07 - - 0s
## 0 0 - 0 - 2.6771e+07 - - 0s
## 0 2 - 0 - 2.6771e+07 - - 0s
##
## Explored 4914 nodes (1716 simplex iterations) in 0.37 seconds
## Thread count was 1 (of 4 available processors)
##
## Solution count 100: 2.67712e+07 2.67732e+07 2.67736e+07 ... 3.14836e+07
##
## Optimal solution found (tolerance 1.50e-01)
## Best objective 2.677119300000e+07, best bound 2.677119300000e+07, gap 0.0000%Sum solutions
marxan_end_ssoln <- marxan_end_soln %>%
mutate(sum = rowSums(.[6:105])) %>%
select(id, cost, status, locked_in, locked_out, sum)
hist(marxan_end_ssoln$sum,
main = "Selection frequencies",
xlab = "Number of runs units were selected")
Create map
Join polygons with output
parcels_marxan_end <- inner_join(parcels, marxan_end_ssoln, by = "id")
ggplot(data = parcels_marxan_end) +
geom_sf(aes(fill = sum),
color = "white",
size = 0.05) +
scale_fill_gradient(low = "slategray2",
high = "navy") +
labs(fill = "Summed \nSolution") +
theme_minimal()
Add basemap
end_spec <-
ggmap(morrobay) +
geom_sf(data = parcels_marxan_end,
aes(fill = sum),
color = "white",
size = 0.1,
alpha = 0.85,
inherit.aes = FALSE) +
coord_sf(crs = st_crs(4326)) +
scale_fill_gradient(low = "slategray2",
high = "navy") +
labs(title = "Endangered & Threatened Species",
fill = "Summed \nSolution",
x = NULL,
y = NULL) +
theme_minimal()
end_spec
Combine maps with patchwork
This step isn’t necessary, but is really convenient if you want to present your maps together in your report. patchwork uses PEMDAS to combine ggplots together into one graphic. For example, plot_1 + plot_2 will make an image of the plots side by side, whereas plot_1 / plot_2 will make an image of plot_1 over plot_2. For more information, see https://github.com/thomasp85/patchwork
First, make the all_spec graph without a legend so that the combined image has only one legend.
all_no_lgnd <- all_spec +
theme(legend.position = "none")Next, combine the graphs with patchwork
spec_graphs <- all_no_lgnd + end_spec
spec_graphs
Finally, save the image to add to your report!
ggsave("spec-graphs.png")